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Abstract Graphical models for finite-dimensional spin glasses and real-world combinatorial optimiza- 
tion and satisfaction problems usually have an abundant number of short loops. The cluster variation 
method and its extension, the region graph method, are theoretical approaches for treating the com- 
plicated short-loop-induced local correlations. For graphical models represented by non-redundant or 
redundant region graphs, approximate free energy landscapes are constructed in this paper through 
the mathematical framework of region graph partition function expansion. Several free energy func- 
tional are obtained, each of which use a set of probability distribution functions or functionals as 
order parameters. These probability distribution function/functionals are required to satisfy the re- 
gion graph belief-propagation equation or the region graph survey-propagation equation to ensure 
vanishing correction contributions of region subgraphs with dangling edges. As a simple application 
of the general theory, we perform region graph belief-propagation simulations on the square-lattice 
ferromagnetic Ising model and the Edwards-Anderson model. Considerable improvements over the 
conventional Bethe-Peierls approximation are achieved. Collective domains of different sizes in the dis- 
ordered and frustrated square lattice are identified by the message-passing procedure. Such collective 
domains and the frustrations among them are responsible for the low-temperature glass- like dynamical 
behaviors of the system. 

Keywords region graph ■ belief propagation ■ Bethe-Peierls free energy • Edwards-Anderson spin 
glass • partition function expansion ■ graphical model 



1 Introduction 

Spin glass is a paradigm for basic research on physics of collective behaviors induced by disorder and/or 
frustration ;2j . It is also a tool box for applied research on hard problems originated from fields outside 
the conventional domain of physics. The mean-field theory of spin glasses was originally developed 
through the replica trick and by assuming the property of replica-symmetry-breaking (RSB) [niH51l5S] . 
It was later re-derived using the cavity method j2E»27^ without introducing replicas. Another advantage 
of the cavity method is its applicability to single problem instances with fixed disorder parameters. The 
cavity method was successfully extended to spin glass models defined on finite-connectivity random 
graphs [25] . confirming and extending the earlier theoretical results obtained by the replica method 
50,29,30]. These theoretical advances paved the way for the fruitful applications of spin glass theory 
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to interdisciplinary problems in computer science, information theory, and biological sciences [13,23 ; 
for example, constraint satisfaction and combinatorial optimization [28,18 , perceptual learning with 
discrete synapses [5], signal transmission and processing [14,47], network structure inference [2"4"li51,3, 
03J, compressive sensing [9TIT7]. 

The most crucial simplification made in the mean-field RSB cavity method is the Bethe-Peierls 
approximation [41138] . Roughly speaking, this approximation ignores all the possible correlations among 
the set of vertices that are nearest-neighbors to a specified central vertex. There have been several 
theoretical attempts to include the effects of correlations to the cavity method 31,37,6,4b. A simple 
and general scheme of partition function expansion was introduced in |52j along this research line, 
following the initial work of Chertkov and Chernyak [3] on Ising spin glasses. By this expansion scheme, 
the mean-field free energy expression at each level of replica-symmetry-breaking can be formally derived 
without imposing any physical assumptions, and the correction contributions to these mean-field free 
energies are expressed as a sum over looped subgraphs. The strong constraint that the subgraphs 
with dangling edges all have completely zero correction contribution to the free energies leads to a 
set of message-passing equations, such as the belief-propagation equation and the survey-propagation 
equation, that must be satisfied by the auxiliary probability distribution functions introduced in the 
expansion scheme. 

For readers not familiar with spin glass theories, the partition function expansion scheme [52] can 
serve as a straightforward mathematical approach to the essentials of the RSB mean field theory. 

Finite-dimensional spin glass models are still very challenging for theoretical studies. One of the 
major reasons is the abundance of short loops and the associated (possibly strong) local correlations. 
Short loops are also presented in numerous graphical models derived from real-world applications 
in optimization, inference, and constraint satisfaction problems. For such spin glass systems, the free 
energy correction contributions from loopy subgraphs may be comparable to the mean-field free energy 
values. As a consequence, the prediction power of the the mean-field theory can be severely damaged. 
A powerful conceptual framework for finite-dimensional lattice models is the cluster variation method, 
which was invented by Kikuchi |15j and later further developed by many authors (reviewed in [341 
[39]). The key insight behind the cluster variation method is that, in a finite-dimensional system, the 
correlation between two vertices decays exponentially with the distance (as long as the system is 
not approaching a continuous phase transition), and therefore short-range correlations dominate the 
statistical property of the system. Short-range correlations among clusters of neighboring vertices 
are considered in a variational way in the cluster variation method [T51IT1I34 1 I39 ] . As an extension of 
the cluster variation method, Yedidia and coworkers proposed a region graph method [53] to tackle 
local correlations in finite-dimensional graphical models with more computational flexibility. A set of 
regions is specified in the region graph representation of a spin glass model [53] . Each region contains 
a subset of vertices and some (or all) of the interactions within these vertices. A partial order can be 
defined among these regions, which is represented by a set of directed edges between pairs of regions. 
A variational free energy is defined on the region graph, the minimization of this free energy leads 
to a set of generalized belief-propagation equations [S3]. Recently, the generalized belief-propagation 
was applied to the two-dimensional (2D) Edwards-Anderson model by Ricci-Tersenghi, Mulet and 
co-workers [81I40U2T] . 

As demonstrated in a very recent paper [56| . the partition function expansion scheme of [52J can 
also be applied to spin glasses in the region graph representation. In the present work we give a 
detailed discussion on how to construct a region graph RSB mean-field theory for spin glass models 
and also to obtain the correction expressions to such a theory. Two sets of message-passing equation, 
namely the region graph belief-propagation equation and region graph survey-propagation equation, 
are derived through the partition function expansion. We also prove that, for non-redundant region 
graphs, the region graph Bethe-Peierls free energy as derived through the partition function expansion 
is equivalent to the variational free energy functional used in |53j . The region graph belief- propagation 
(rgBP) equation is applied to the 2D square-lattice Ising model and the 2D square-lattice Edwards- 
Anderson (EA) spin glass model, and our numerical results are compared with the results obtained by 
the conventional belief-propagation (BP) procedure. 

The prediction on the Ising model's Curie temperature by the rgBP is better than the predic- 
tion of the conventional BP. The exactly known transition temperature can be approached if region 
graphs with sufficiently large maximal regions are used in the rgBP. When applying the rgBP to the 
square-lattice EA model at sufficiently low temperature, a most interesting numerical observation for 
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us is that, the rgBP identifies many small collective domains of vertices in the square lattice (see 
figure If). Vertices in each of these collective domains are strongly correlated (and probably change 
states collectively), while the remaining vertices outside of all these small localized domains serve as 
paramagnetic background. Such highly heterogeneous patterns are in some respect similar to the pat- 
terns of dynamical heterogeneity in dense liquid when approaching the glass transition (reviewed in 
10, 1 2]). A key difference is the collective domains in the square lattice do not move in space, only their 
areas increase as temperature decreases. These collective domains bring in many new time scales to the 
system's low-temperature dynamics, they should be essential to understand many of the fascinating 
low-temperature glass-like behaviors of this system. 

For the square-lattice EA model, the existence of many localized collective domains above the 
paramagnetic background naturally points to a way of improving the power of the rgBP. The basic idea 
is to use larger maximal regions for the identified collective domains, while the sizes of the maximal 
regions in the remaining parts of the lattice keep to be as small as possible. This hybrid adaptive 
strategy will be explored in a future work. We hope it will be helpful in achieving deeper insights on 
the heterogeneous (but globally still paramagnetic) behaviors of this and other 2D spin glass systems. 

The next section introduces the general graphical model and its factor graph and region graph 
representations, and briefly reviews the cluster variation method. Section [3] derived the region graph 
replica-symmetric (RS) mean-field theory as well as the region graph belief-propagation equation; the 
connection with the generalized belief-propagation of [53] is also discussed here. Section [1] derived 
the region graph first-replica-symmetry-breaking (fRSB) mean-field theory and the associated region 
graph survey-propagation equation. The RS mean-field theory (rgBP) is then applied to the 2D Ising 
model and the 2D Edwards- Anderson model in Sec. [5j We conclude this work and list some possible 
future projects in Sec. [6] The paper have four appendices. 



2 The Model and Graphical Representations 

Consider a very general model system of N particles (z = 1,2,..., N) interacting with each other and 
with the environment. Each particle i is fixed in space and therefore has no translational degrees of 
freedom, it is fully characterized by an internal state variable Xi. For example, in Ising models, the 
state variable is binary, xi = ±1; in Potts models, the internal state can choose among Q different 
values, Xi € {1, 2, . . . , Q}; in Heisenberg models, Xi is a three-dimensional continuous vector of unit 
length. In this paper we assume, for notational simplicity, that the state variables Xi are discrete and 
can take only a finite number of different values. The microscopic configuration (state) x of the whole 
system is defined by the states of all its particles, x = {x\, X2, ■ ■ ■ ,xn}- The configuration energy H(x) 
is supposed to have the following additive form 

N M 

H{x) = Y i E i {x i ) + Y i E a {x da ). (1) 

i—1 a—1 

On the right side of ([I]), the first term is contributed by external forces (fields), with each energy 
depending only on a single particle i (if i is free of external forces, then Ei(xi) = 0). The second term 
is contributed by M internal interactions (a = 1, 2, ... , M), each having an energy E a . Let us denote 
by da the set of particles that are involved in interaction a and by x_Q a = {xi\i € da} a joint state 
of the particles in this set. The internal interaction energy E a is a function only of Xg a . For example, 
in a two-body Ising model an internal energy has the form E a = —J a XiXj with J a being a coupling 
constant, then da = {i,j} and x da — {xi,Xj}. 

After equilibrium is reached in an environment of temperature T, the probability of the system 
being in a configuration x obeys the Boltzmann distribution 

where (3 = l/(fcgT) is the inverse temperature, ks being the Boltzmann constant. The partition 
function Z((3) is expressed as 

N M 

Z(P) = £ e -**fe) = J2H^i) II ^feea)> ( 3 ) 

x x i—1 a—1 
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where ipi and ip a are, respectively, the Boltzmann factor for the external field on particle i and the 
internal interaction a, 

fcfc) = e -^ x <\ ^fejj = e-f* E -<3»»\ (4) 
From Z(j3) we can define the equilibrium free energy F(/3) as 

F{0) = -ilnZCS) = Y J p B{x)H{x) + k B TY J PB{x)\nP B {x). (5) 

X X 

The expression ([5]) for the equilibrium free energy can be extended to a general non-equilibrium 
situation. Given an arbitrary probability distribution p(x) on the configurations of the model ([I]), the 
Shannon entropy functional is defined as [7] 

S[p] = -k B Pfe) Wfe), ( 6 ) 

X 

and the free energy functional F[p] is 

F[p] = J2p(^H(x)-TS[p]. (7) 

X 

It is easy to check that the absolute minimum of F[p] is equal to the equilibrium free energy F(f3) and 
this minimum is achieved only when p(x) = P B (x) . 



2.1 The Factor Graph Representation 

The model defined by Q can be represented conveniently by a factor graph G of variable nodes 
(shown as circles, denoting the N particles k, . . .), function nodes (shown as squares, denoting the 
M internal interactions a, b, c, . . .), and edges between pairs of variable and function nodes (i, a)F]The 
factor graph G is a bipartite graph: all the edges are between a variable node and a function node, 
and an edge (i,a) exists between a variable node i and a function node a if and only if variable i is 
involved in the internal interaction represented by a. A comprehensive review on factor graphs can be 
found in |19) . 

As a simple illustration, we show in Fig.[T]part of the factor graph for the two-dimensional Edwards- 
Anderson spin glass model [TT] on a square lattice. The energy function of the EA model is defined 
as 

where the internal energies are due to ferromagnetic (Jy > 0) or anti-ferromagnetic (J^ < 0) spin 
couplings along the edges of the square lattice, and h® is the external magnetic field on particle i. 
We use (Ji G {— 1, +1} instead of Xi to denote the binary spin state of particle i. The spin configuration 
of the whole system is a_ = {<j\, <j-x, . . . , cat}. 



2.2 A Brief Summary of The Cluster Variation Method 

The cluster variation method (CVM), as originally proposed by Kikuchi [15], is an approximate method 
to calculate the Shannon entropy ^k. It decomposes S[p] into contributions from different clusters of 
variable nodes (see [T] for an easy-to-access introduction to CVM) . A cluster in the CVM is defined as 
a non-empty subset of the N variable nodes. The total number of possible clusters is 2 N — 1. For two 
clusters C\ and C2, we say that G\ < C2 (and C2 > C\) if and only if C\ is a subset of C2 (in cases 



1 We follow the convention in the literature and use letters k,l, . . . to denote variable nodes and letters 
a, b,c,d, . . . to denote function nodes. 
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Fig. 1 The factor graph for the two-dimensional Edwards- Anderson model. Each variable node % (circle) 
represents a lattice site i, it has a spin state <7i and an energy Ei{cri) = —h1<7i caused by an external field 
ft?; each function node a (square) between two variable nodes i and j represents a spin coupling with energy 

E a {tJi,(Tj) = —JijOiOj. 



that C\ < C 2 and also C\ > C2, then C\ = C2). If C\ is a strict subset of C2, we denote as C\ < C 2 
(and C2 > Ci). This partial ordering corresponds to the following function £(Ci,C2) 

C(Cl ' C2 ) = { 0, othe^wife. 6,2 ' (9) 
This function can be regarded as a matrix, and it is invertible: 

C(Ci, C 2 ) M (C 2 , Os) = XI A*(Ci, <? 2 )C(<? 2 , C 3 ) = 5(Ci, C 3 ), (10) 

where <5(Ci,C2) is the Kronecker delta such that <J(Ci,C 2 ) = 1 if Ci = C 2 and = otherwise. The 
function n(C\,C2) is the Mobius inversion function [15] defined as 

v ' [0, otherwise. 

In the above equation, |C| denotes the number of elements in cluster C. 

The configuration of a cluster C is defined as the collection of states of all its elements and is 
denoted as x c = {xi : i £ C}. Given a probability distribution p(x) for all the TV variable nodes, the 
marginal distribution for the configuration x c is 

PcisHc) = J2 Pfe)- (12) 

Following ([6]), we can define for each cluster C its Shannon entropy functional as 

S C [pc] = -kB^2pc(x c )\np c (x c ). (13) 

Let us define an entropy increment ASc for each cluster C as (see, for example, pQ) 

AS C = ]T (-l) |Ch|C '^c =^S C >H(C,C)- (14) 

C'<C c 



Then from (14) and (10 1 we obtain that 

S c = J2 AS c<(C',C) = Y As c- (15) 



C'<C 
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In other words, the entropy Sc of each cluster C is the sum of the entropy increments of all its 
sub-clusters (including C itself). 



If we take C as the set of all the variable nodes in G, the above equation (151 then states that, 
the entropy functional S[p] of the whole system is the sum of the entropy increments of all the 2^ — 1 
different clusters, i.e., S[p] = ASclpc]- There are an exponential terms in this summation. The 
insight of Kikuchi [15] was that, the entropy increment ASc often decays quickly with the size \C\ of 
the clusters, therefore a good approximation of S[p] can be obtained by summing over only a subset 

of small clusters. Let us choose a set of maximal clusters c[ m \ C% , . . . , C p m \ These p clusters are 
maximal in the sense that any of them is not a sub-cluster of any another cluster in this chosen set. 
These p maximal clusters and all their sub-clusters form a set, denoted as P. We then have the following 
approximate expression for the entropy functional: 

S[ P ] « V AS c [pc]. (16) 



CEP 



This is the CVM approximation. 



Inserting ( 14 1 into the above expression, we get 

S[p] « a cSc[pc], (17) 

CEP 

where the coefficient ac of a cluster C is expressed as 

oo=5>(C,C)= £ (-1)I C 'H C I (18) 

C'EP {C'\C'EP,C>C} 

It is easy to check that each of the p maximal clusters (say C^) has the coefficient a C ( m ) = 1. For 
other clusters of P, their coefficients can be calculated iteratively as 

a c = l- ac '- ( 19 ) 

{C'\C'EP,C>C} 

The applications of the CVM method were reviewed in [54T39] . This method can be combined with 
Suzuki's coherent anomaly method |46| to compute the critical exponents of a given finite-dimensional 
system. 



2.3 The Region graph Representation 



In the CVM approximation (16), for each maximal cluster, all its sub-clusters are included into the 
cluster set P. This requirement was relaxed in the work of Yedidia and co-authors, who proposed a 
more flexible region graph representation [53] . A region graph R is formed by region^ and a set of 
directed edges between the regions. A region a of a factor graph G includes a set of variable nodes and 
a set of function nodes, with the condition that, if a function node a belongs to a region a, then all 
the variable nodes in the set da also belong to the region a. A configuration of a region a is defined 
by the states of all the variable nodes in this region, x a = {xi\i £ a}. 

If there is a directed edge pointing from a region p to v in the region graph R, then it must be 
the case that v is a sub-region of p (that is, all the variable nodes and function nodes of region v also 
belong to region p). We use p, — > v to emphasize the directness of the edge; sometimes if the direction 
of an edge is unimportant for the discussion, the notion (p, v) is also used. It should be emphasized 
that, a region v being a sub-region of another region p does not necessarily indicate the existence of 
a directed edge from p to v. If there is a directed edge from region p to v, then p is regarded as a 
parent of v and v a child of p. If there is a directed path from region a to region v, we say that a is 
an ancestor of v and v is a descendant of a. The ancestor-descendant relationship between two regions 
a and v is denoted by v < a and a > v. (If v is not a descendant of a, such a relationship does not 
exist even if v is a sub-region of a.) 



2 In this paper, we use Greek symbols a, f3, 7, . . . to denote the regions of a region graph R. 
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Fig. 2 A non-redundant region graph R for the factor graph shown in Fig. ^ Only a part of the full region 
graph is shown. There are three types of regions: each 'square' region (e.g., region a) contains n x n variable 
nodes and 2n(n — 1) function nodes (in this example, n = 2), its counting number is c = 1; each 'rod' region 
contains n variable nodes and n — 1 function nodes, its counting number is c = — 1; and each 'stripe' region 
(e.g., region fj) contains n x 2 variable nodes and n function nodes, with counting number c = 1. Each stripe 
region is connected to two rod regions, and each rod region connects to one stripe region and one square region. 
The short red lines between two regions indicate the parent-child relationship (the directions of these edges 
are not shown, as they are obvious). In this region graph, variable node i appears in 5 different regions, and 
function node a appears in 3 different regions. 



Each region 7 is assigned a counting number c 7 , which is determined recursively by 

C-y 1 ^ ^ Ca • 

{a|a£i?,a>7} 



(20) 



Notice the similarity between (20) and (19). In the region graph R, the subgraph Ri induced by a 



variable node i is defined as the subgraph that include all the regions containing i and all the directed 
edges between these regions. Similarly, R a denotes the region graph induced by function node a, it 
includes all the regions containing a and all the directed edges between these regions. The region 
graph R and its associated counting numbers {c 7 } are required to satisfy the following region graph 
conditions [53] : 

(1) For any variable node i, the induced subgraph Ri is connected, and the sum of counting numbers 
within this subgraph is unity: 

E c 7 = L ( 21 ) 

76-R; 

(2) For any function node a, the induced subgraph R a is connected, and the sum of counting numbers 
within R a is unity: 

E c i = L ( 22 ) 

Given a factor graph G, one may construct many different region graphs R, all of them satisfying 
the region graph conditions. A region graph R is considered as being non-redundant if it satisfies the 
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following additional condition: the region subgraph Ri induced by any variable node i contains no 
loops (it is a tree of regions) . Region graphs which do not satisfy this tree condition are referred to as 
being redundant. As we will discuss in detail, non-redundant region graphs have some nice properties 
for partition function expansion. 

A simple region graph R is shown in Fig. [2] for the two-dimensional EA model. Although this region 
graph contains many loops at the region level, it is non-redundant as can be easily checked. 



3 Partition Function Expansion on a Region Graph 

Given a region graph R for a factor graph G, we now expand the equilibrium partition function Q 
as a sum over the contributions of region subgraphs of R. The mathematical framework of |52] will be 
followed. 



3.1 Introducing auxiliary states of variable nodes and removing redundancy 



First of all, because of the constraints (21 1 and (22 1, the partition function can be expressed as 



x aeR 



(23) 



The microscopic configuration of each region a is denoted as x a = {xf \i £ a}, where xf is the state 
of variable node i in region a. For a function node a in region a, the variable states at its neighborhood 
{xf \i £ da}. A Boltzmann factor <P a (x a ) is defined for region a as 



are denoted 



^(£a) = n^«)II^(^a)- 



(24) 



aea 



A variable node i may belong to two or more regions (say a, 7, . . .). If this is the case, we regard the 
states of node i in the different regions, xf,x],..., as different variables. With the introduction of 
these new variables, the partition function expression (|23|) is re-written as 



{x a } a£R 



N 



{x a } aGR 



■yeRi 



(25) 



The Kronecker delta functions 8(xf,Xi) guarantee that, the partition function is contributed only by 
those region configurations {a; Q |o! € R} for which the states xf,x],... of each variable node i in 
different regions are equal to the same value. 

Let us for the moment focus on the region subgraph Ri induced by variable node i. Assume Ri 
contains n$ regions. The expression within the square brackets in (251 requires that the states of 
variable node i take the same value among the rii regions, which is equivalent to (rii — 1) constraints. 
These (rii — 1) constraints of each variable node i can be conveniently implemented in the following 
way: (1) Attach on each edge (/1, v) of region graph 7? a variable node set, the cross-linker set 
which is initially empty. (2) For each variable node i € {1,2,..., N}, choose (n, — 1) different edges of 
the induced region graph Ri in such a way that the rii regions of Ri form a connected tree with these 
edges; then for each of these (rii — 1) chosen edges (/z, v), add i to its cross-linker set \i#v. (3) After all 
the cross-linker sets are constructed, for each variable node j in any cross-linker set , a constraint 
5(&j, crj) is applied, forcing the states of j to be the same in the two connected regions ju and v. 

The region graph shown in figure [2] is non-redundant, and there is no need to remove redundancy. 
The left panel of figure [3] shows a redundant region graph R for the 2D Edwards- Anderson model. This 
region graph has three types of regions, and each variable node (such as node 1) appears in 9 different 
regions. The region subgraph R t induced by variable node i contains 12 directed edges. One particular 
way of removing redundancy is indicated by the X and XX symbols on the directed edges of R. After 
removing redundancy, the resulting graph is shown in the right panel of figure [3] This graph 
is actually equivalent to the region graph shown in figure [2] (we will return to this point later). 
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Fig. 3 Removing redundancy from a region graph. (Left panel) a redundant region graph R for the factor 
graph shown in Fig. [TJ There are three types of regions in this region graph, as indicated by the three different 
colors. Each variable node (e.g., i,j,k,l) appears in 9 different regions. A X symbol on a directed edge fi — » v 
(between a parent region fi and a child region v) means that the single common variable node (say i) of fx and 
v is not included into the cross-linker set (that is, = 0). A XX symbol on the edge /j, — > v means 

that the two common variable nodes (say i, j) of /j, and v are not included into the set fJ,#v (again, fi#v = 0). 
(Right panel) the resulting non-redundant graph i?* after all the redundancy in R has been removed. 



By this construction, it is easy to check that the states of each variable node i is associated with 
(n-i — 1) constraints. For some edges (/x, v) of the region graph R, the cross- linker set might be empty, 
fi#v = 0. We can remove all such edges from the region graph R, and the resulting graph is denoted 
as R# (see the right panel of figure p5l for an example) . The partition function ( 25 ) is re-expressed as 



e n 

{x a } aGR 



n 



n * 

iEfj,#u 



(26) 



We should emphasize that the above-mentioned procedure of removing redundancy does not change 
the counting numbers, it only assigns a cross-linker set \i$v (might be empty) on each directed edge 
(/i, v) of the region graph R. 



3.2 Introducing auxiliary cavity probability functions 



For each edge (//, v) of the region graph RF we now introduce two auxiliary probability distribution 
functions, p^u (x"n u ) and Pv^^(^t^t v )- They are non-negative and normalized but otherwise arbi- 
trary. The microscopic configuration jJJ^, denotes the states of variable nodes of set [i^v in region fi. 



i.e., x^jtv = { x i N ^ A*#^}; similarly, x£# u = S ^#^1- We further rewrite (26) as 

Z[fi) 



e n ^fe.) n Pi- 

{x a } a L 7£3#a 



n 



P^v (ZuftfJPv^-H (^i#v) 



where d#a denotes the set of nearest neighboring regions of region a in region graph RF . 
It is helpful to define two partition function factors as follows: 



X 



E 



(27) 



(28) 
(29) 
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Then (27) is written as the following simple form 



where the coefficient Zq is equal to 



and -4(^,1/) is expressed as 



z =nz a n ~7~~' 



a£R {p,v)£R# 



y^—tv \±v#n IPv^yp \±n#u I 

and w a (a; Q ) is a probability distribution defined as 



(30) 



(31) 



(32) 



(33) 



We regard ^(^j,) as small quantities and expand the edge-product of (30), and finally obtain 

i+ y, a 



Z(/3) = Z 



rCR# 



(34) 



where r denotes any subgraph of the region graph R& , which contains a subset of the edges of R& and 
the associated regions. The correction contribution L r of a region subgraph r is expressed as 



Lr= Yl w i^y) II A (^M^#v^l#n)- 



(35) 



3.3 Region graph belief-propagation equation 

Many of the subgraphs of the region graph R& contain dangling edges. A dangling edge (/i, v) in a 
subgraph r is such an edge that if it is cut off from subgraph r, either region [i or region v (or both) 
will be an isolated region of r. For a dangling edge (/x, ^) of subgraph r, suppose the region /x has only 
this edge attached to it. The correction contribution L r of subgraph r is calculated to be 



Pfi^-u&vftft) Ylix^ #v ^>x" Tli£fi#v $( x i> x i)P^p( x ^ v )Pn-^v( x u^ l J 
.Pti-^vi^v^ii) ^2tx* Ex^ #fi Ylien^v 3( X i ' X i)P^f( X -fj,#u)P^^(^-u#n) 

where p^v{g%,4t v ) is a probability distribution function determined by 



E^ug n ?w^# T ) n 



(36) 



iW(^ #i J = S^,(K^|7 G ^W) (37) 

(38) 



In the above expressions, d& n\v denotes the set of nearest-neighboring regions of ix in the region graph 
J?*, but with v being removed from this set. 
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From (36) we know that, if the auxiliary probability distributions {p^-H/iPu-t/j,} are chosen as a 
the following equation, 



fixed-point o 



n iv^GeM n 



(39) 



then all the subgraphs r of the region graph i?# with at least one dangling edges will have zero 
correction contribution. Then the free energy of the system is expressed as 



F{P) = F -~hx 



I 1 > lJ T loa 

r lo °PCR# 



(40) 



where r loop denotes a subgraph of R# that is free of dangling edges (in r loop , each region has at 
least two edges attached). Equation ( [39] ) is referred to as the region graph belief-propagation (rgBP) 
equation. It ensures that all the subgraphs of with one or more dangling edges have zero correction 
contribution to the partition function Z(/3), therefore greatly reduces the number of terms in the free 
energy correction expression. 



It is not difficult to check that, the expression of the rgBP equation (39) for the region graph RF 
shown in the right panel of figure [3] is equivalent to that of the rgBP equation for the non- redundant 
region graph shown in figure [5J In this sense, i? # of figure [3] (right panel) is equivalent to R of figure [2J 



3.4 Region graph Bethe-Peierls free energy F 



When all the loop correction contributions in (40) are neglected, the remaining term F is an approx- 



imation to the true free energy F(f3). We refer to Fq as the region graph Bethe-Peierls free energy, its 
expression is 

F = £/a- £ fM> ( 41 ) 

with 



£*«go n 

P"/—HX\2La#j) 



(42) 



The free energy F can also be regarded as a functional of the probability functions {p^utPv^,^}- A 
nice property of this functional is that, the partial derivative of Fq with respective to any probability 
function p^ v vanishes at a fixed point of the rgBP equation (39). This can easily be checked by 
showing that 



E 2V 



(43) 



where the probability function p u -¥u has been defined through (38). 



From the variational property of Fq, we know that each fixed point of the rgBP equation ( 39 1 locates 
a stationary point of the functional Fq and vice versa, namely there is a one-to-one correspondence 
between the stationary points of the Fq functional and the fixed points of the rgBP equation. A 
stationary point of F could be a minimum, it could also be a maximum, or be a saddle point. 

Using the set of probability functions {p^ v ,p v ^^) as order parameters, the functional Fq gives 
an approximate description of the free energy landscape of model 0. If Fq has only one station- 
ary point (expected to be a minimum), then the rgBP equation has only one fixed point (solution). 
This simple situation can be referred to as the 'replica-symmetric' (RS) case, in connection to the 
physics approaches based on the replica method and the cavity method [271I23] . In many non-trivial 
problems, however, the Fq functional may turn out to be very rugged in shape and have many sta- 
tionary points, then (39) has multiple solutions. This later complex situation can be referred to as the 
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Fig. 4 The region graph Bethe-Peierls free energy Fo as a functional of probability functions may have only 
a single stationary point (left panel), or it may have multiple stationary points including minima, maxima, 
and saddles (right panel). This qualitative change in functional shape is induced by variation in environmental 
temperature or other control parameters. 



'replica-symmetry-broken' (RSB) case. We will discuss this RSB case in the next section|4j A schematic 
illustration of the qualitative change in shape of Fo is shown in Fig. [4j 

The free energy Fq can be expressed as a sum over the contributions of all the regions. For any 



region a of the region graph R, we have the identity c a + J2 1>a c 7 



a 

^ CaF a 



7>a 



a>fi 



1. Then Fq can be re- written as 

(44) 



Each region a contributes a term c a F a to Fq. The expression of F a is 

-^a " Jot 



/(".TO + 



(45) 



Notice that (45) has the same expression as (41), the only difference is that the summations are 
restricted to the subgraph containing region a and all its descendants. Therefore, F a can be regarded 
as the region graph Bethe-Peierls free energy of this subgraph a. 



3.5 Factor graph as the simplest region graph 

It has been demonstrated in |52j that the conventional belief-propagation (BP) equation can be derived 
from the framework of partition function expansion. Actually BP is just a limiting case of the more 
general rgBP equation. 

Given a factor graph G, let us consider the following simplest region graph R: there are N 'variable' 
regions denoted by i = 1, 2, . . . , N, each of which contains a single variable node i, with counting number 
Ci = 1 — ki (ki = \di\ being the degree of node i in G); and there are M 'function' regions denoted by 
a = 1, 2, ... , M, each of which contains a function node a and the set da of all the nearest-neighboring 
variable nodes of node a, with counting number c a = 1. Such a simple region graph is non-redundant, 
although it in general contains many loops. 

Let us define two new probability functions bi^ a (xi) and b a ^i(xi) as 

^(x^p^ajXj) jj-\ Xi )p a ^{Xi) 



Using the rgBP equation ( [39] ) , one can easily check that the newly defined probability functions satisfy 
the conventional BP equation: 

A(xi) n 

) V'afegJ II b j->a(Xj) 
h f n _ c€di\a XaaVi jGda\i 

x c£di\a x 0a j<£da\i 
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The probability (a;,) can be interpreted as the state distribution of variable node i in the 
absence its interaction with function node a; while the probability b a ^i{xi) can be interpreted as the 
state distribution of i if it only interacts with function node a. 

The free energy F for this simplest region graph reduces to the conventional Bethe-Peierls free 
energy. Its expression is F = J^aeG ^ a — EieG^* ~ 1)-^- The free energies Fi (for a variable node i) 
and F a (for a function node a) are expressed as 



F a = 



hi 



(48) 



Only the messages b a ^i(xi) from the function regions (parents) to the variable regions (children) 
appear in the above two expressions. We now proceed to demonstrate that, such a property holds for 
a general non-redundant region graph. 



3.6 Parent-to-child message-passing in a non-redundant region graph 



In a non-redundant region graph i?, the subgraph i?, induced by any variable node i is a connected 
tree. Because of this property, i? # is identical to R and the cross-linker set for each directed edge fj, v 
contains all the variable nodes in v (/i#^ = £ v}). It is easy to see that, if there is a directed path 
pointing from a region a to another region 7, then such a path must be unique in R; furthermore, the 
subgraph R a induced by any function node a is also a connected tree. 

When we approximate the free energy F{0) by the region graph Bethe-Peierls free energy Fo, the 
corresponding approximate expression for the marginal configuration distribution of a region 7 is 



Pa- 



•7(^7#a)' 



(49) 



This marginal distribution has the following consistency property: if (/1, v) is an edge of the region 
graph i?, then 

Y Pm(^m) = Y P»(2v)> ( 50 ) 
where x^ u = {xi\i £ (J-ftv} is the configuration of the variables in the cross- linker set /z#za Equation 



( 50 1 ensures that the marginal configuration distribution of the cross-linker set \i^v is the same whether 
it is inferred from (x„ ) or from p v (x v ) . 

When R is non-redundant, we have d^j — d^y in ( 49 ) , where d"f is just the set of nearest-neighboring 



regions of region 7. Then, with the help of (39), the expression (49) is rewritten as 



(51) 



In this expression, the region set J 7 contains 7 and all its descendants, i.e., 2~ = {77I77 < 7}; B 1 is 
another region set with the property that any region in B 7 does not belong to 7 7 but is parental to at 
least one region in 7 7 , see figure [5j The configuration x v for any region v £ I 1 should be understood 
as x v = {x]\i E v). 

For a non-redundant region graph R, we can prove (see the appendix [A]) that, the marginal prob- 
ability distribution (|51[) is equivalent to 



-da 



a£7 



iG7 



n m A 

{(/i->-i/)|/iGB 7 ,vel~,} 



(52) 



where rn^ v {x_ l/ ) is a probability distribution that satisfies the self-consistent equation (53). The mes- 



sage m^^fij is interpreted as the (cavity) probability that the variable set {i £ v\ takes configuration 
x u in region /x when all the interactions that also appear in region v are not considered (namely, setting 
^i(xi) = 1 for any i £ v and ipa(x da ) = 1 for any a £ v, while all the other interactions of region 
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Fig. 5 For any region 7 in a non-redundant region graph R, the set 7 7 is formed by 7 and all its descendant 
regions, and the set _B 7 is formed by all the other regions that do not belong to J 7 but are pointing to regions 
of J 7 by one or more edges. Each white disk in this figure represents a region. 



H are not modified) . The marginal probability expression ( 52 1 has the nice property that only the 
parent-to-child messages m M _>.„ are needed, but not the messages from children to parents. 



The self-consistent equation for mfj,^. v (x v ), derived from (39) after a lengthy process (see Appendix 
fAI), is intuitively simple: 



n 



(53) 



For a non-redundant region graph R, the marginal probability distributions Pj{x^) should satisfy 
the following edge-consistency condition, namely 



V(/i -> 1/) e iJ. 



(54) 



To check this consistency, plug ( 52 1 into the two sides of the above equation, we get the follow equation 
that must be satisfied by m^ u : 



n n i>a(xa a ) n u^) 



V-^fe,) — C- 



n to 

{(a-^laes^e/^vr,,} 



a— ^-77 V— 77 



{■n^x\rie{B v ni l j.)\iJ,,xei v \u} 



(55) 



where the normalization constant C is determined by m ll ^ u {x y ) = 1. For a non-redundant graph 

R, the set (-B„ n Ip,)\n is empty, and hence the denominator of the expression on the right side of (55 1 
is equal to 1. Since (|55|) is equivalent to (53 1, the edge-consistency condition (54 1 holds. 



If the region graph R is redundant, the marginal probability distribution of each region 7 can 
no longer be written in the form of ( |52[ ) , and the edge-consistency condition ( 54 ) is no longer valid 
(but the weaker condition (50) is still valid). In [5"3"|, however, Yedidia and co-authors defined that 



a region 7 has a marginal probability distribution ( |52[ ), and then they required that this marginal 
probability distribution should satisfy the edge-consistency condition (54), which leads to the self- 



consistent equation (55) for m^lij. Such a treatment might be able to get good results for some 



problems, but in our opinion it lacks a solid theoretical fou ndation. 

For a non-redundant graph R, the free energy F a (45) can also be expressed by parent-to-child 



messages (see Appendix [B] for details). The final expression is 
1 



F n 







In 



n 



v(x_ v ) 



(56) 
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The free energy expression ( 56 ) also appeared in [53) as a basic assumption (both for redundant 
and non-redundant region graphs). A theoretical challenge for redundant region graphs, worth to be 
explored in future work, is to derive from the framework of partition function expansion an expression 
for F a using only parent-to-child messages. 



4 Region Graph Expansion for the Grand Partition Function 



If the region graph Bcthc-Peierls free energy functional (41 ) has only a single stationary point, a global 



minimum, then the rgBP equation ( 39 1 has only a single solution (fixed point) . Reaching this unique 
fixed point is computationally not hard. One may perform message-passing iteration along the edges 
of the region graph R, or one may work on the functional Fq and find its unique minimum by gradient 
descend methods. The situation is a little bit more complicated, for example, if the Fq functional has 
a single minimum and one or more saddle points. In this situation, we expect that the rgBP equation 
still has only one stable fixed point, which can be reached by direct iteration (with some damping to 
accelerate convergence) or by minimizing Fq. 

The real non-trivial situation occurs when the Fq functional has multiple minimal points (see the 
right panel of figure [4]) , and consequently the rgBP equation ( 39 1 has multiple stable fixed points. We 



take Fq as an approximate free energy landscape for the system under study. Then each minimal point 
of Fq can be regarded as a metastable state s of the system, it describes approximately the system's 
collective behavior within certain timescale t s . The timescale r s is determined by the shape of the free 
energy landscape surrounding metastable state s, especially the relative height of the closest free energy 
saddle point. For some mean-field models (such as spin glasses on random graphs [25] ) r s is expected 
to approach infinity in the thermodynamic limit N — > oo (see, e.g., 32,22|), then a metastable state 
s is a true thermodynamic Gibbs state (i.e., ergodicity is broken at N — > oo). For finite-dimensional 
systems, the timescale t s of a metastable state s may be finite even in the thermodynamic limit. 



In this section, we regard each fixed point of the rgBP equation ( 39 ) as a macrostate, no matter 
whether it is a minimum of Fq or a maximum or saddle point. Viewing each macrostate as an energy 
level, it is natural to define a grand partition function E as pMf2l)lll5"5"Il5"2"] 

S(y;0) = J2 exp[-2/^ (s) (/?)], (57) 

s 

where .Fq (/3) is the value of Fq(/3) at the rgBP fixed point s; the parameter y is the inverse temperature 
at the level of macrostatesj^The grand partition function E(y; /3) is associated with a grand free energy 
G(y; f3) through 

At a given inverse temperature y, the statistical weight of a macrostate s is 

% 

Eexp[-yF (s) (/?)] 



G(y;p) = --lnS(y;p). (58) 
V 



exp[-yF (a) (/?)] 

V b(s) = — - r 777737' ( 59 ) 



From Vb(s) the mean value of free energy is expressed as (Fq(/3)) v = J2 S '^ > b{s)Fq {0), Arid the 
entropy density £{y; (3) of macrostates is expressed as 

E(y^) = -~Y l V B (a)hiVB(s). (60) 

S 

£{y\ ft) is also called the complexity in the spin glass literature [25]. It is easy to verify that G(y; j3) = 
(F (P)) y -fS{y;P). 

3 Although (|57j) sums over all the stationary points of Fo, we expect that at sufficiently large values of y the 
grand partition function will be dominantly contributed by the free energy minimal points. 
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We are interested in obtaining an approximate expression for the grand free energy G(y;j3). For 
this purpose we follow again the framework of partition function expansion |52j . First notice that 



n /d?^ /d? 5{Pn^v - B^ y )5{p v ^^ - B v ^)t 



-yF 



(61) 



where J Dp^^ means integrating over all possible probability measures Pfi-^vix^ ) on the edge (/i, v) 
of the region graph R#. Because of the Dirac function 8{p^ y — B li ^, u {{p 1 ^, ll \^ G d^fj,\i/})), only the 
fixed points of ( 39 ) contribute to the grand partition function. 

We then introduce on each edge (fj,, v) of the region graph RW two probability measures Pu-y v {pu,-^v) 
and Pv^^pv^p). Using the expression (41 ) for the free energy Fq, the grand partition function is re- 
written as 



s( y; /3)=n n 



5{p^u - B^ v )8{p v . 



.(62) 



This expression is very similar in form to (27 1, therefore the method of partition function expansion 
can be directly applied to ( 62 1 . As a result we obtain that 

1 



G(y;P) = G --hx 

y 



i+ E 

r i°°PCR# 



The first term on the right of ( 63 ) is expressed as 



where 



9-i 



= In 



9{^, V ) 



Go(y; P) = E f7 - E 



Y[ J Dp a _, 7 P a _,. 7 (p a _, 7 )e yf -> 



In 



Dp^vDpv^P^vip^^P^^pv^^e a/ <^> 



(63) 
(64) 

(65a) 
(65b) 



The grand free energy Go as expressed by ( 64 ) can also be regarded as a functional of the probability 
functional {P^ v {p^ v ), P^^p^^)}. 

The correction contribution Lr of a subgraph to the grand free energy has a similar expression 
as (35) [52] • To ensure that any subgraph with at least one dangling edge has vanishing correction 
contribution to the grand free energy, each probability function Pu-^-v (Pu->^) needs to satisfy the 
following equation 



3i£ ,.\ ,. V / 



where 



n 



■yed#fi\ 



v J~Dp 1 ^ fl P 1 ^ tl {p. 



7~W 



-Vfu 



I fl—H> 



P 



In 



n 



(66) 



(67) 



Equation ( |66[ ) is called the region graph survey-propagation equation, in correspondence to the 
1RSB mean- field theory of spin glasses 25,28J. If we neglect the loop correction contributions to 
G(y;f3), then the grand free energy functional Go gives an approximate description of the system's 
free energy landscape at the level of macrostates. It can again be verified that, the first variation of 
Go with respect to any of its arguments P M _).„(p M _>.„) is identically zero at a fixed point of (66). 



17 



If we neglected all the loop correction contributions to G(y;0) and approximate it with Go, then 
an approximate expression for the mean free energy of macrostates is 



<W)>, 



d(yG 
dy 



E 



n 



-yfy 



n 



J Dp. 



E 



'V—tfl-Pfl—tV {Pfl—t- 



u)Pv->»(pu-*)f(n,v)e- vf M 



And the complexity E(y; /3) is expressed as 



y((F (f3))y - Go) 



N 



(69) 

j3, and to determine the maximal 



It is particularly interesting to determine the complexity value at y 
value of y at which the complexity becomes negative [33J . 

It is not a simple task to solve numerically the region graph survey-propagation equation ( 66 ) . A 
direct approach is to iterate ( 66 1 on the region graph, with each probability functional rep- 
resented by a sample set of probability functions P/j,-^u(^,^u)- A ma j° r complication is the rewcighting 
factor e~ v ^^" in (66). In the special situation of y = (3, the reweighting can be replaced by introducing 
new auxiliary probability functions |22U18j . Several reweighting tricks were discussed in the thesis of 
Zdeborova |54| . which may be helpful for solving the region-graph survey-propagation equation for 
general y values. 



The region-graph survey- propagation equation ( 66 1 applies to a general redundant or non- redundant 
region graph R. Rizzo and co-authors |40j derived a set of generalized survey-propagation equations 
using the replica cluster variation method, which are different from the region-graph survey-propagation 
equation (66). The connection between these two approaches needs to be further studied. 



5 Numerical Results on the Two-Dimensional Ising and Edwards- Anderson Models 

We now apply the region graph belief-propagation equation to the 2D Ising and Edwards-Anderson 
models, mainly for the purpose of testing its performance. Both models are described by the energy 
function For the Ising model, all edge coupling constants Jij are equal to a positive value J; for 
the EA model, are independent and identically distributed random variables, taking value J and 
— J with equal probability. All the external fields h® in Q are set to be zero for simplicity. In the 
numerical calculations, J and the Boltzmann constant ks are both set to be unity, so the energy unit 
is J and the temperature unit is J/fcs. 

We consider L x L square lattices with L coupling interactions in each dimension. If we assume 
periodic boundary condition on both dimensions, the total number of vertices (variable nodes) in the 
lattice is N = L 2 ; if open boundary condition is used on both directions, the total number of vertices 
isA=(L + l) 2 . 

The region graph for the square lattice is constructed to be non-redundant. It has three types of 
regions, the 'square' regions, the 'stripe' regions, and the 'rod' regions. Each square region contains 
n x n vertices, each stripe region contains n x 2 vertices, and each rod region contains nxl vertices. 
All the coupling interactions between the vertices of a given region are also included into this region. 
Each square region is connected to four rod regions, each stripe region is connected to two rod regions, 
and each rod region is connected to a square region and a stripe region. The region graph R at n — 2 
is shown in figure [2j where a square region and a stripe region have the same shape. For the general 
case of n > 2, if we regard each square region as a 'giant vertex' and each stripe region plus its two 
connected rod regions as a 'giant bond', then the region graph is again a square lattice of giant vertices 
and giant bonds. In this sense, our region graph representation is a coarse-graining of the original 
lattice that keeps its topology unchanged. 

The partition function of the 2D Edwards- Anderson can be calculated exactly by polynomial algo- 
rithms (see, for example, Ref. [48] and references therein). Although the heuristic rgBP message-passing 
approach can only obtain an approximate value for the free energy, it is very efficient in estimating all 
the N single-variable marginal probabilities simultaneously. 



18 




©-Eh-© 

Fig. 6 A local part of the region graph R shown in figure^ The counting numbers for a square region a, a 
stripe region fi, and a rod region a are, respectively, c a = 1, = 1, c a = —1. 



5.1 Region graph belief-propagation equations at n = 2 

A subgraph of the region graph R at n = 2 is plotted in figure [6] For notational simplicity we denote 
the square regions and the stripe regions by Greek symbols (such as a and ^i) and denote a rod region 
just by the index of its function node (such as a and c). Consider a square region a and a rod region 
a in figure [6] The two probability distributions between these regions can be parameterized as 

Pa-+a{<n,<rj) « exp(ph^ a a % + f3h^l a a 3 + /3J^l} a a % a 3 ), (70a) 
Pa-+a{<ri,<rj) oc exp(/?^ Q a l +/3^ a (T J +/3J^ Q a l( 7,). (70b) 

Similarly the two probability distributions between a stripe region \x and a rod region a are expressed 
as 

p^ a (ai,aj) oc exp^/iW^a-j + /3h^ a aj + PJf*l} a <7iO-j), (71a) 



The self-consistent equations for these set of parameters are derived from the rgBP equation |39|) and 
listed in appendix [C] 

A t rivia l solution of the rgBP equation is the paramagnetic one with the fields in the expressions 
( 70a H 71b) all being identically zero, 

ft&.=*&.=0, fc& a = fc& a = 0, fc<& =*<& =0, ^=4^=0. (72) 

The stability of the paramagnetic solution is analyzed through a set of linearized rgBP iterative equa- 
tions listed in appendix [d] All the fields such as ha^-a arL d h^U a are randomly initialized, and their 
values then evolve according to the linearized rgBP equations. If all the fields finally decay to zero, 
the paramagnetic solution is stable, otherwise it is unstable and the original rgBP equation has other 
stable fixed points. 

The number (2™ — 1) of parameters needed to completely characterize a probability distribution 
of the rgBP equation grows quickly with the number n of vertices on a boundary line of the square 
region. The paramagnetic solution and its stability analysis have been worked out up to n = 10 for 
the Ising model. For the EA model we have only considered the simplest case of n = 2. 



5.2 The ferromagnetic Ising model 



At sufficiently high temperatures T the paramagnetic solution ( 72 ) is the only solution for the rgBP 
equation at n = 2. This trivial solution becomes unstable at the critical temperature j'j" -2 ) — 2.65635 
(periodic boundary conditions). This value is higher than the exact transition temperature T c = 2.26919 
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Fig. 7 The critical temperature T at which the paramagnetic fixed-point of the rgBP equation becomes 
unstable (periodic boundary conditions). The integer n is the number of vertices on a boundary line of a 
square region. The exactly known phase-transition temperature T c = 2/ln(l + y/2) ~ 2.2692 is indicated by 
the horizontal dashed line. The red solid curve is a fitting function T = 2.2376 + 0.6875n~ ' 7140 to the data 
with n > 2. 



|16j . but a little bit lower than the value of Tc — 2.88539 as obtained through the conventional 
belief-propagation equation (i.e., n = 1) [Hl3"B"]. 

The paramagnet-ferroniagnet transition temperature as predicted by the rgBP equations decreases 
if larger square regions are used. Figure [7] demonstrates that the predicted critical temperature Tj™^ 
can be fitted by the following curve 

T^=T™ + - h , (73) 



with T c °° = 2.2376±0.0037, a = 0.6875±0.00f 7, and b = 0.7140±0.0f 08. The fitted value T c °° is slightly 

lower than the exact critical point T c — 2.2692. The data can also be fitted well by Tc = T c + a 1 /n b , 
with a' = 0.688 ± 0.007 and b' = 0.81 82 ± 0.0078. 

The free energy density, mean energy density, entropy density, and magnetization as a function of 
temperature T are shown in figure [8] and compared with the exactly known results of Onsager |35j . 
As n increases, the results are closer to the exact value. The free energy as obtained by the rgBP 
equations is an upper bound to the true free energy value of the system, but the difference is small at 

(n) 

n > 2. Both the energy density and the entropy density have a kink at the rgBP critical point T c , 
but this kink becomes more and more weaker as n increases (in the limit of n —¥ 00 the exact Onsager 
solution should be reached). Using the theoretical framework of coherent anomaly method the 
results obtained at different values of n can be used to predict the critical exponents of the 2D Ising 
model. The results of such an exercise (to be carried out) will be reported elsewhere. 

Under the periodic boundary condition, the instability temperature of the paramagnetic solution 
shown in figure [7] is independent of lattice side length L. But this is not the case for the open boundary 
condition. Under the open boundary condition, the shorter the lattice side length L, the more stable 
the paramagnetic solution is (see figure [9]). However, this difference in threshold temperature between 
the periodic and the open boundary conditions becomes very small for L > 100. The rgBP equation 
at n = 2 have two stable ferromagnetic fixed points even for a very small 5x5 square lattice (open 
boundary conditions) if the temperature is lower than 1.92 (see figure[9]). For larger open square lattices, 
the ferromagnetic solutions are stable at higher temperatures. Of course for finite lattices there is no 
real phase transition. The two low-temperature ferromagnetic solutions are interpreted as describing 
the two metastable states of the square lattice. They are stable because the correlation length in the 
system exceeds the length scale of the maximal square region of the region graph. If n exceeds the 
correlation length of the system, the paramagnetic fixed point of the rgBP equation will again be 
stable. 
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Fig. 8 Results obtained on the ferromagnetic Ising model (periodic boundary conditions) using the conven- 
tional BP and the rgBP at various values of n. The results of Onsager's exact solution are also shown for 
comparison, (upper left) free energy density; (upper right) mean energy density; (lower left) entropy density; 
(lower right) magnetization. 
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Fig. 9 Threshold instability temperature of the paramagnetic solution for the Ising model on a L x L square 
lattice with open boundary condition. Circular points are obtained by the BP approximation, square points 
by the rgBP at n = 2. At each value of L, the threshold temperature as predicted by rgBP at n — 2 is lower 
than that predicted by the conventional BP. 
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5.3 The Edwards- Anderson model 

A set of single instances of the 2D EA model with ±J coupling constants are randomly generated. 
For each instance the threshold stability temperature Tc f the paramagnetic solution of the rgBP 
equation at n = 2 is determined numerically, and this value is compared with the threshold value 
as obtained by the conventional BP approximation. As shown in figure |10[ the performance of rgBP 
(n = 2) is better than that of BP, but Tc H j s still positive, and it increases with the side length 
L of the periodic square lattice. Figure [10] also suggests that the finite-size corrections to the critical 
temperature Tc decrease as 1/ InL. This finite-size scaling behavior is different from the observation 
of 0. 

The EA model on a square lattice has no real spin glass phase at finite temperatures (see, for 
example, 33 J44j). At any positive temperature T the system has only a paramagnetic phase, and the 



magnetization on each vertex is equal to zero in the long time limit. Apparently, the results of figure 10 
with the paramagnetic solution being unstable at positive temperatures are in contradiction with the 
absence of finite-temperature spin glass phase. 

This apparent discrepancy can actually be removed. The instability of the paramagnetic rgBP 
solution does not mean the system is in a spin glass phase. It just signifies the emergence of some 
collective domains in the square lattice, the length scales of these collective domains exceed the length 
scale (= n) of the region graph's maximal square region. The detailed arguments go as follows. 

An elementary square (including four coupling interactions) of the square lattice is called a plaque- 
tte. There are two types of plaquettes, frustrated or non-frustrated. A plaquette is said to be frustrated 
(non-frustrated) if the product of the four edge coupling constants on its boundary is negative (posi- 
tive) [15] . It is an obvious fact that the four edge coupling energies of a frustrated plaquette can not be 
simultaneously minimized. For the EA model studied in this paper, on average one-half of the plaque- 
ttes in each problem instance are frustrated, and these frustrated plaquettes are randomly distributed 



on the 2D lattice (see the upper left panel of figure 11 for a concrete example) 



The abundance of frustrated plaquettes destroys the long-range ferromagnetic correlations in the 
system. However, the local density of frustrated plaquettes fluctuates considerably at different parts 
of the lattice. There are some small patches of the lattice that are mainly formed by non-frustrated 



plaquettes. For example, in the 64 x 64 periodic square lattice of figure 11 (upper left panel), there 
is a 5 x 5 patch centered at position (62, 28) in which only 4 of its 25 plaquettes are frustrated. The 
spin coupling interactions in such small patches are essentially ferromagnetic in nature (under a gauge 
transformation of the spin variables and the coupling constants [49]). Some of these ferromagnetic 
small patches may considerably exceed the maximal n x n square region in size, and at low enough 
temperatures, the ferromagnetic correlation length will exceed n, making the rgBP fixed point to be 
locally ferromagnetic (after the gauge transform) in these patches but paramagnetic in the remaining 
parts of the square lattice. Given a patch with a specified contour shape and a specified density of 
unfrustratcd plaquettes, the probability to discover such a patch is higher in a square lattic e w ith 
longer side length L. Therefore, the observation that yj" -2 ) increases with lattice size L (figure 10 ) is 
consistent with figure [9j 

The marginal probability distribution of the two spins <7j and o~j of a rod region a (see figure |6| is 
calculated as 

Pa{0-i,<Jj) OC e f3Jtj<7i ' 7i p a ^ t a{0-i 1 0-j)p l _ l ^ t a{0-i,0'j) 

(x exp[/3(^ a + h% a )ai + + + + " ^K^] ■ (74) 

Then we get the magnetization of vertex i as 

_ tanh[/3(fei^ a + h%g)] + tanh[/3(fo^q + fe<& 8 )] tanh[/3(jff&, + jj$ a - Jg)] 
mi ~ 1 + tanh[/3(/^A> a + h ( ;h a )] tanh[0(/i&„ + J#2>«)] tanh[/3(J^ a + J<j$„ - J -)] ' 



(75) 



The magnetizations of all the vertices in the square lattice can be obtained in a similar way. We define 
the mean absolute magnetization of a plaquette as the average of the absolute magnetizations of its 
four vertices. 



Figure 11 shows the pattern of mean absolute plaquette magnetizations at a given temperature T, 



for a single instances of the EA model on a 64 x 64 periodic square lattice. For f) < 0.60177 the rgBP 
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Fig. 10 The threshold instability temperature of the paramagnetic solution of the conventional BP and the 
the rgBP equation at n = 2 for L x L periodic square lattices. Each data point is obtained by averaging 
over 10 single instances of the Edwards- Anderson model. The dashed lines are two fitting curves of the form 



T 



bj log 2 L, the two fitting parameters a, b are shown in the figure. 



equation has only the paramagnetic solution, and all the plaquette mean magnetizations are zero. At 
(3 = 0.606, the mean absolute magnetizations of a small domain of the square lattice are nonzero 



(figure 11, upper right panel). As temperature further decreases, this collective domain enlarges in 
area, and also other collective domains start to form (figure [IT] lower right panel, (3 — 0.65). As 
temperature further decreases, d iffer ent collective domains start to get into contact and they compete 
for the boundary vertices (figure 11 lower left panel, j3 = 0.75)Q 



The heterogeneous patterns and its evolution shown in figure |11| i n some respect are similar to 
the phenomenon of dynamical heterogeneity in structural glasses [T0lll2j . This link deserves to be 
explored more deeply. The existence of many collective domains and the frustration effects between 
these domains very probably are responsible for the glassy-like low-temperature dynamics of the 2D 
Edwards- Anderson model. 

Figure [TT] also suggests a way to improve the rgBP prediction power. We infer that different local 
regions of the square lattice have different correlation lengths. To consider more precisely the corre- 
lations in the collective domains, a conceptually easy way is to construct region graphs with larger 
regions for these collective domains, while small regions are used for the remaining parts to lower 
computational complexity. This adaptive strategy needs to be implemented in future work. 

We end this section by emphasizing that, the BP and the rgBP iterative process are still able to 
converge when the paramagnetic fixed point becomes unstable. For the 2D square-lattice Edwards- 
Anderson model (periodic boundary conditions), we found that BP converges as long as T > 1.52, in 
agreement with earlier simulation results [2"U1I5]. The rgBP at n — 2 converges at even lower tempera- 
tures. 



6 Conclusion and outlook 

In this paper we gave a detailed description of the region graph partition function expansion approach 
(first introduced in |56|). and obtained approximate expresses for the free energy and generalized free 
energies of a general graphical system with an abundant number of short loops. A series of message- 
passing equations (such as region graph belief-propagation and region graph survey-propagation) were 
derived in the expansion process. We have applied the rgBP equation to the square lattice Ising model 



4 At such low temperatures, the rgBP equation starts to be difficult to converge, probably because of the 
frustration effects due to domain competitions. 
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Fig. 11 Emergence of domains with collective dynamics, (upper left) Distribution of frustrated plaquettes 
in a single instance of the 2D Edwards-Anderson model with periodic boundary conditions. Each frustrated 
plaquette is shown as black, while each non-frustrated plaquette is shown as white. For this instance, the 

paramagnetic solution of the rgBP equation at n = 2 becomes unstable at yj™- 2 ) — 1.6618 (/3 = 0.60177). 

(upper right) At temperature T = 1.6502 (f3 — 0.606), which is just slightly below yj™- 2 ^ one small collective 
domain begins to emerge, as demonstrated by non-zero magnetizations for some plaquettes in a connected 
cluster of plaquettes. The color of each plaquette encodes the mean value of absolute magnetizations as averaged 
over its four vertices, (lower right) At a further decreased temperature T = 1.5385 (/3 = 0.65), another collective 
domain becomes quite evident. As the temperature further decreases, more collective domains form, and the 
formed collective domains enlarge in size and their boundaries start to be in contact. As an example, the mean 
absolute magnetization of each plaquette is shown at T = 1.3333 (/3 — 0.75) at the lower left panel. 



and Edwards- Anderson model and found that it outperforms the conventional belief-propagation equa- 
tion. We also demonstrated that the fixed points of the rgBP equation reveal the heterogeneous pattern 
of collective domains in the square lattice. An adaptive strategy of improving the rgBP performance 
was also suggested. 

As discussed in the previous sections, there are many issues remain to be explored. The adaptive 
rgBP scheme for the 2D Edwards- Anderson model needs to be implemented. Another very interesting 
issue is the effect of adding redundancy to the region-graph. A redundant region-graph have been used 
for the 2D Ising model and Edwards- Anderson model in various papers within the framework of cluster 
variation method 15,34.53,8,21]. The performance of rgBP needs to be tested on such a redundant 
region-graph. We have performed some preliminary computations and found that the results of rgBP 
depend on the particular way of removing redundancy. For example, the region graph R# shown in 
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the right panel of figure [3] actually is equivalent to the non-redundant region graph R of figure [2j Some 
other ways of removing redundancy are able to make the paramagnetic rgBP fixed point to be stable 
at even lower temperatures. 

For a general redundant region graph R, the rgBP equation ( 39 1 is not equivalent to the generalized 



belief-propagation equation (55 1 of [53) . We are working on the issue of deriving the generalized belief- 
propagation equation from the approach of partition function expansion. 

The three-dimensional Edwards-Anderson model is believed to have a true spin glass phase. This 
system will be studied using rgBP and the region-graph survey-propagation equation in a future work. 

Particles in a dense liquid have translational degrees of freedom. It remains to be seen whether a 
similar partition function expansion scheme can be worked out for systems with mobile particles. If 
approximate free energy landscapes can also be built for such systems with the help of message-passing 
equations, it should be very helpful for understanding structural glasses and supercooled liquids. 

A different theoretical approach (the replica cluster variation method) of studying finite-dimensional 
spin glasses has been explored in }4"01[51l!?l") . The message-passing process of this replica cluster variation 
method appears to be much more complicated than the simple rgBP process. 



A Derivation of parent-to-child message-passing equation 



Here we give a detailed derivation of equations (52 I and (53 1. These equations are valid for a non-redundant 
region graph R. We assume 7? to be non-redundant in this whole section. 

First, the region subgraph formed by region 7 and all its des cend ants i s a connected tree (see the blue- 
shaded area of figure This property ensures the equivalence of (51 1 with (491. Notice that a region fi £ 73 7 
may point to two or more r egio ns of the set 7 7 



Applying the definition (241 and then the two identities (211 and (221, we obtain that 



n *w = n[^fej]" si ' ,nRa "nhM* 61 "^ * 



(76) 



The set Rb\I-, contains regions of subgraph Rb except those also belonging to set J 7 , and similarly for Rj\I-y 
Because Rb and Rj are two connected tree subgraphs, we have 



Y c "= Y Y Y c ° 

r,6-R b \/ 7 veR b nI-y {(/i->!/)|M6B 7 } aGR^" 



Y c "= Y Y Y c ° 



(77) 



In the above equation, R%~*" denotes the branch of the tree Rb that is still connected with /i if the directed 
edge fi — > v is removed; and Rt^ v has the same definition, i.e., it is the branch of the tree Rj that contains 
region but not v. The possibility that a region /i £ 73 7 might point to two ore more regions in 7 7 does not 
affect the validity of (77 1. The reason is simple: if v\ and V2 are two children of /1 in J 7 , then v\ and vi do not 
share any fun ction node nor any variable node in common. 

Based on ( |77[ ), ( |76| l and ( |51| l, we obtain the important expression ( |52[ ). In that equation, the parent-to- 
message m IJ ,^{x v ) is defined as 



£ 



(78) 



6Gi^ j£v 

up to a normalization constant (to be fixed by £^ ?7l^-> !/(£„) = 1). 

Using the expression (391 for the probability distribution p fl ^. u (x u ), it is easy to show that 



E 



X 



y n 



n p a ^\{x x ). 



(79) 
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Notice that 



n = n Im^jV^ 1 ^ n kMl^^ 
n k] 1 °e"°v» Ca n kl 1 °™< c " n kl""'^^-^ n k-l" eRin<w " )C ' 7 - 



(80) 



'-LE[i\u 



Then we have 



E 



n i> a n v> 



nW aeHr " nkl ° 6fl ^" n 
n kl - 6 ^ c ° n kl q ^ 
n kl ° e ^ u " c ° n k- 



E 



E c Q + E 

qefl n(J„\J„) 



E 



It is easy to check that, for a function node c 6 1/ £ 



(81) 



(82) 



where R'^ v \I ll denotes the set formed by all the regions in the region subtree R^" except those which are 
also members of the region set 7 M . Similarly, for a variable node k £ v £ 1^, we have 



(83) 



with R^ fI/ \I fl being the the set formed by all the regions in t he region subtree R£ " except those which are 
also members of the region set I M . With these two equalities, (81 1 is re- written as 



E 

nkl ° eRr " nkl aeR > " n 



E c Q 



n i>« n in 

.a£fi\v i(EfJ,\u 



x n k]^ e ^ c °nkl ^^ Ca n kl ae ^ c< *nkl ^ R ^ I - C \ (84) 



Combining the above equation with ( 78 1 leads to 



n 



E 



E 



n kl • ,6 ** v " nkl * 6R ^"^ n nkl" e * 

6G/A" c6f L J {(a^A)|aeS^,Agi" M \/„} dSA L J 



E c v 



n kl nk 



E 



E c >7 



{(«— >-A) |a£ 



n nkl" 6 ^ 

B^AG/fA-Tn} iSA L J 



(85) 



Using the the properties ( 77 1 for the region trees Rb (induced by each function node b G n) and Rk (induced 
by each variable node k 6 fx ), it is not difficult to verify that the expression in t he curly brackets of the above 
equation is equal to 1. Therefore we arrive at the message-passing equation ( 53 ( for m^-^ v (xj)- 



B Derivation of the free energy expression (56) 



We demonstrate that, for a non-redund ant region graph R, the region graph free energy Fq can be expressed 
as F = Ecgjj c a F a , with F a given by ( 56 1 . 

First, we notice that /( M ,„) as defined by (|42| can be expressed as 



/W) = fu + g In 



(86) 
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where fv-yp is defin ed through ( |67[ ). In writing down this equation, we have assumed that fi is a parent of v. 
From equation ( 45 1 we then obtain that 



Fa — fa + ^ ^ fv- 



(87) 



On the other hand, based on the definition ( 42 1 for f a we derive that 



e n n p^-fe 

£a l6 f o {(/i->")lMe-B a ,y6-T a } 



E /- 

{(/^-^JlMe/a} 



From the last two expressions we then get the following simple formula for F a 



F a = -^ln 



e n n p/^fej 

E„ 17 €■!"<* {(p->i/)|//£B a , l/6-fc,} 



(89) 



This formula is very similar to (561, but not yet identical 



Now we replace Pn-tvixJ) of 1 89 1 by rn^-n/(sD u ) through the relation (78 1, and obtain that 



e -HEnhM'^^nM 3 *)] n&Rr ^p^ 



(90) 



We need to prove that 



E 



E 



e^ e ^ En KM ,,e " r " nKM " 6 " r " 



= 0. (91) 



To prove this, it is first noticed that the left side of (91 1 is equivalent to 



{(/J->i/)} v {a\ii£B a ,i/eI a } 



E E c >, 



e e ^ta^nKM "* r " nKM 



j£v 



(92) 



For a non-redundant region graph i?, the sum in the curly brackets of the above expression is identical to zero, 
i.e., for each directed edge fj, —+ v. 



E Ca ~ E c " ~ E °^ ~ 1 ~ 1 ~ °" 



(93) 



Combining equations (44 1, (89 1 and (91 1, we obtain the objective equation 



F o = J2 Ca {~l ln 



(94) 
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C Self-consistent equations at n — 2 

For the local structure shown in fisure [6] the square-to-rod messages between the regions a and a are: 



= h 



,0) 



j; 



(0 



+ 



hi 



cosh[^(fti + J d )] 



2/9 Lcosh[/3(ft ; - J d )] 



+ 



■hi 



1 + tanh[/3 J c ] tanh[/3(/ii + Jd)] tanh[/3(/i fc + J 6 )] 



4/3 L 1 + tanh[/3 J c ] tanh[/3(ft i - J d )] tanh[/3(ft* - J b )} 



1 + tanh[/3J c ] tanh[/3(hi + J d )] tanh[/3(/i* - J b ) 



1 + tanh[/3J c ] tanh[/3(ft; - J d )] tanh[/3(/i fc + J&) 



= h 



0) 



+ 



2/3 



In 



cosh[/3(fe fc + J,)] 
cosh[/3(/i fe — Jb)] 



tanh[/3 J c ] tanh[/3(/i fc 



1 + tanh[^ J c ] tanh[/3(/i fc 



+ - 



■hi 



1 + tanh[/3J c ] tanh[/3(h fc + J,)] tanh[/3(ft; - J d ) 



— Jii + 
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In 



4/3 Ll + tanh[/3 J c ] tanh[/3(/i fc - J 6 )] tanh[/3(ft; + J d ) 
1 + tanh[/3 J c ] tanh[/3(fe; + J d )] tanh[/3(fe fc + J b )] 
1 + tanh[/3 J c ] tanh[/3(^ + J d )] tanh[/3(fe fe - J b )] 

1 + tanh[/3J c ] tanh[/3(h; - J d )] tanh[/3(/ifc - J,) 



.1 + tanh[/3J c ] tanh[/3(/i; - J d )\ tanh[/3(/i fc + J 6 ) 
where we have introduced several shorthand notations 



+ j b )} tanh[/3(ft; + J d )] 



J 6 )]tanh[/3(/ii - J d )] 



(95a) 
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Similarly, the stripe-to-rod messages between the regions fi and a are: 
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(■inn) 
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(96a) 



(96b) 



(96c 



On the same edges (a, a) and (/i, a), the rod-to-square and rod-to-stripe messages are much simpler and 
are expressed as 



•,(*) 



= ft' 1 ) 
= h 
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(0 
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D Stability analysis of the paramagnetic solution at n — 2 



At the paramagnetic fixed point (72 I, the effective coupiings such as Ja^la and are determined self- 

consistentiy through (95c I and (95c I. Then the rgBP iteration equations for the fieids are iinearized. The 
coefficients of the iinearized equations are obtained by the foliowing expressions: 



!■(«) 



dh™ f _ tanh[/3J im ] (l - tanh 2 [/3 J™ a ] tanh 2 [/3J jn ]) 
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9ft« . 
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(98c 
(98d 
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(98f 



(98 
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The iinearized rgBP equations for the fields are iterated on a given region graph. During each sweep of 
the iteration, the output field messages of each square region of the region graph are updated once, and the 
maximum among the absolute values of all the updated fields is recorded. If this maximum decays to zero with 
the iteration sweeps, the paramagnetic fixed point is then declared as stable. When the paramagnetic solution 
is unstable, this maximum value will eventually increase with iteration sweeps (after a transient decreasing 
stage). 
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